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A detailed description is provided of a numerical model of the microcavity optical parametric 
oscillator (OPO). The model is solved on a two dimensional grid, in the time domain, using an 
alternating-direction implicit method. Results are presented for a system pumped with a gaussian 
spatial profile beam. It is shown that at high pump powers the OPO signal forms a ring, with low 
emission at the centre where the pumping is strongest. 
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I. INTRODUCTION 

The polariton physics in a strongly-coupled semiconductor microcavity includes interesting non-linear optical effects. 
When resonantly pumped with an off-axis laser, optical parametric oscillator (OPO) behaviour can be obtained, with 
the spontaneous appearence of signal and idler emission^. Since the excitonic non-linearity is x > the phase matching 
conditon requires 2ui p = u) s + tOi and 2k p = k s + fcj, which means that the pump, signal and idler can all be close to 
resonance with the polariton dispersion. This behaviour has been investigated theoretically by a number of authors, 
using analytio&a and numerical^ models. However, all these models differ from the experimental setup by assuming 
that the pump is a plane wave, with uniform excitation of the system, and uniform polariton fields. While such 
treatments provide useful insights into the OPO, a fuller understanding of the observed behaviour requires a theory 
which takes into account the finite excitation spot used in experiments. The purpose of this paper is to describe a 
numerical treatment of a two dimensional microcavity polariton system pumped with a Gaussian beam, solved in the 
time domain on a discrete spatial grid. 

II. THE NUMERICAL MODEL 

The model for the OPO is based on the non-linear optics treatment of Refi It consists of coupled two-dimensional 
cavity and exciton fields, (j) c and X , with a local x non-linearity in the exciton to represent the effects of short-range 
interaction. Both fields have dispersions, described by effective masses M c and M x , and are damped with widths 70 
and 7x- The coupling is characterised by Q, the Rabi splitting, and A, the detuning at normal incidence. The cavity 
field is driven by a spatially dependent external pump field f p (r) exp (iujpt). The time evolution of the fields is thus 
determined by 
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To be exact, the form for the x^ nonlinearity should be k</> x . However, for a numerical treatment, it is better 
to chose re|</> x | 2 </> x , as this is independent of the origin of energy, which can thus be shifted to the exciton frequency. 
The modified form correctly gives the scattering effects of interest here, but does not describe frequency multiplying 
effects, which would require extremely small time-steps in the numerics. 

For the numerical treatment, this model is discretized on a two dimensional grid, with the differential operators 
are replaced by hopping terms. The solution is obtained in the time domain, as an initial value problem, using the 
Crank-Nicholson method^. This leads to a semi-implicit time-step formula of the form 

(1 + \iH8t) = (1 - \iHSt) , (3) 

where the vector ip^ = (<j> c , <fi x ) T describes the cavity photon and exciton amplitudes on the site (ij) at time-step n. 
The precise form of the matrix operator H will be specified below. 

The problem with solving Eq.JSJ) in two dimensions is that H is a large matrix (size iV 2 x N 2 , where N is the length 
of the edge of the grid), whose structure does not lend itself to efficient numerical solution: there are non-zero terms 
far from the diagonals. It is better to use an operator-splitting treatment, the alternating direction implicit method^, 
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writing H = H x + H y , where H x 
divided into two sub-steps: first, 



is solved for every j, then 



and H y respectively contain only x and y hopping terms. Then, each time-step is 
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for every i. Each of these one dimensional strip problems is independent, so it is necessary to solve only 2N (N x N) 
sets of equations. Moreover the matrices H x and H y contain non-zero terms only on the diagonal and first three 
sub-diagonals, which permits a very efficient solution by factorisation^. 
The part of H x involving coupling to site i is 
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where each H is a 2 x 2 matrix 
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The hopping terms are t^; = (2Af c ^) _1 and t£ = (2Af x (S^) _1 , with S x the ^-direction grid spacing. The value of </> x 
used in the non-linear term is the old value at that site, ie from the right hand side of Eq.QJ or Eq.(|SJl. H y has a 
similar form, except 5 y replaces S x , though in practice an isotropic grid is used. The factors of ^ multiplying 7J slte 
are needed to avoid over-counting, because they occur in both H x and H y . 

Eq.QJ, and similarly Eq.JSJl, is most efficently solved in two stages: first an intermediate variable ip* 1S calculated 
using 

\ (1 + \iH x 5t) = (8) 

then it is straightforward to show that i(j n+ i =tp t — ip n . 

To incorporate periodic boundary conditions, Eq.© has to be modified by adding hopping terms in the off-diagonal 
corners of the matrix, coupling the sites at either end of a strip: 
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Adding these terms breaks the useful diagonal form of the matrix, so a direct solution of Eq.@ would be inefficient. 
However, as there are just four extra non-zero entries, the problem can be handled at little numerical cost using the 
Woodbury formula^ to incorporate the extra terms. is written as 



HZ = HL 



E 

fe=l,2 



(10) 



where the Uk and Vk are chosen to generate the end hopping terms. An efficent choice is 
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where —0i and — @2 are the first two elements along the 
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main diagonal of H x . H' x is a slightly modified version of 
H x : the u(£)v terms produce, as well as the required hopping terms, four additional elements along the main diagonal 
which have to be subtracted from H x to obtain the correct . 

Evaluating the Woodbury formula requires additional solutions of Eq.(JHJ with u\ and u 2 on the right hand side, but 
this is very inexpensive once the factorisation of H x is obtained. Hence adding periodic boundary conditions makes 
little difference to the numerical costs. 
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FIG. 1: Grey-scale plots of the OPO states for Gaussian pump beam at 16°with various powers, (a) and (e) are the signal and 
idler intensities for a peak power of 0.2, just above threshold. The following rows represent pump powers of 0.3, 0.4 and 0.5, 
with the signal states on the left, the idlers on the right. The lengths of the edges of the squares are 51.2/mi, and the FWHM 
of the pump beam is 17.05/im. 

III. RESULTS 

The results described in this section are for a microcavity with physical parameters typical of a good GaAs structure: 
the Rabi splitting, = 5 meV, and the line-widths of the exciton and cavity mode are 7x = 7c = 0.25 meV. The 
structure is on-resonance at normal incidence (k = 0), and pumped at 16° with a Gaussian beam whose FWHM is 
17.05/im. The cavity mode effective mass is taken to be 3 x 10~ 5 m e , and the exciton mass 0.5m e . The spatial grid 
consists of (2 7 x 2 7 ) points with a separation of 0.4^m, and the time step is 12 fs. 

The calculation is started with the fields set to zero, and the solution is then propagated until a steady state is 
reached, typically after about 10 4 time steps or ~ 100 ps. The pump, signal and idler are then separated by storing 
the data for a series of time steps, Fourier transforming to the frequency domain, and applying band-pass filters. This 
process can be carried out in a continuous way using a rolling time window, which allows the observation of dynamical 
effects occuring on time-scales slow compared to the inverse of the mode separation. 

The calculated images show the spatial dependence of the steady-state cavity field intensities in the signal and idler 
modes at various pump powers. It should be noted that they represent the internal cavity field, which cannot be 
directly observed experimentally. The emitted light will undergo some additional spatial filtering due to the small 
/c-width of the polariton branches. 

Fig.l shows the calculated signal and idler images for the steady state at various pump powers above threshold. 
Close to threshold, (a) and (e), only a small spot near to the centre of the image is bright. This is easily understood: 
for a Gaussian beam profile, the threshold will be reached first where the pump is strongest. At higher pump powers, 
the emission spot expands, but, instead of the centre becoming brighter, the spot turns into a ring with low intensity 
in the middle. It can be explained qualitatively using an analytic plane-wave model where the contributions of the 
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signal and idler to the polariton blue-shifts are taken into account*^. The shifts cause the OPO to switch off above 
a certain power, which depends on the incidence angle and detuning of the pump beam. This happens first at the 
centre of the spot, leading to the dark middle. In a Gaussian beam, there will always be a region where the pump 
power is the range to allow the OPO, so an expanding emission ring should be expected as the power is increased 
above threshold. 

Of course, the strong spatial dispersion of the polariton means that a such a simple picture, relating the OPO 
intensity to the local value of the pump, is not completely accurate, but it seems to give a reasonable guide. The 
dispersion effects are responsible for the fact that the ring is not perfectly symmetrical, and breaks up significantly 
at higher pump powers. However, a bilateral symmetry is always mainatincd, defined by the incidence plane, which 
cuts the figure along the a;-axis. 

Finally, at higher powers than those shown in the Fig.l, the numerics predict that a steady state may not be 
reached. Instead, the signal and idler undergo a periodic oscillation, maintaining the general ring shape, but with the 
maximum intensity moving from side to side. 

IV. SUMMARY 

It has been shown that a full understanding of the microcavity OPO requires a two-dimensional treatment of the 
spatial structure of the pump beam and the polariton fields. A numerical model describing the spatial and temporal 
evolution of the system has been derived, and used to predict that the emission spot should form a ring at high pump 
powers. These results demonstrate that experimental measurements without spatial resolution provide only a limited 
picture of the OPO behaviour. 
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